global detritus fish_1 fish_2 P phyt zoopl
global m_phyt m_zoo m_f1 m_f2 m dec pr pr2 in_P uptake_P out_P photo outfl gr;

detritus = zeros(366,1);
detritus(1) = 4;

fish_1 = zeros(366,1);
fish_1(1) = 5;

fish_2 = zeros(366,1);
fish_2(1) = 0.02;

P = zeros(366,1);
P(1) = 0.1;

phyt = zeros(366,1);
phyt(1) = 8;
 
zoopl = zeros(366,1);
zoopl(1) = 0.15;



m_phyt = zeros(366,1);
m_zoo = zeros(366,1);
m_f1 = zeros(366,1);
m_f2 = zeros(366,1);
m = zeros(366,1);
dec = zeros(366,1);
pr = zeros(366,1);
pr2 = zeros(366,1);
in_P = zeros(366,1);
uptake_P = zeros(366,1);
out_P = zeros(366,1);
photo = zeros(366,1);
outfl = zeros(366,1);
gr = zeros(366,1);
   
 



dt = 1;
qv = 0.0125;

for t = 1:366
    if t~=1
        detritus(t) = detritus(t-dt) + (m_phyt(t-dt) + m_zoo(t-dt) + m_f1(t-dt) + m_f2(t-dt) - m(t-dt) - dec(t-dt)) * dt;
        fish_1(t) = fish_1(t-dt) + (pr(t-dt) - pr2(t-dt) - m_f1(t-dt)) * dt;
        fish_2(t) = fish_2(t-dt) + (pr2(t-dt) - m_f2(t-dt)) * dt;
        P(t) = P(t-dt) + (in_P(t-dt) + m(t-dt) - uptake_P(t-dt) - out_P(t-dt)) * dt;
        phyt(t) = phyt(t-dt) + (uptake_P(t-dt) + photo(t-dt) - outfl(t-dt) - gr(t-dt) - m_phyt(t-dt)) * dt;
        zoopl(t) = zoopl(t-dt) + (gr(t-dt) - pr(t-dt) - m_zoo(t-dt)) * dt;
    end
    
    gr(t) = zoopl(t) * 0.4 * (phyt(t) - 0.1) / (phyt(t) + 0.5);
    pr(t) = 0.125 * fish_1(t) * (zoopl(t)-0.1) / (zoopl(t)+0.3);
    pr2(t) = 0.1 * (fish_1(t) - 0.1) * fish_2(t) / ( (phyt(t)+detritus(t))*(fish_1(t)+0.1) );
    
    
    m_phyt(t) = 0.1 * phyt(t) + 0.35 * gr(t);
    m_zoo(t) = 0.1 * zoopl(t) + 0.35 * pr(t);
    m_f1(t) = 0.015 * fish_1(t) + 0.35 * pr2(t);
    m_f2(t) = 0.01 * fish_2(t);
    m(t) = 0.004 * detritus(t);
    dec(t) = 100 * m(t) + qv * detritus(t);
    in_P(t) = qv * 0.8;

    uptake_P(t) = 0.05 * phyt(t) * P(t) / (P(t) + 0.01);
    out_P(t) = qv * P(t);
   
    photo(t) = uptake_P(t) * 100;
    outfl(t) = phyt(t) * qv + 0.05 * phyt(t);

end